icr_parser::icr_parser(string infilename, string outfilename, string dataname){
   infile = infilename;
   outfile = outfilename;
   datafile = dataname;
   cout<<"ICR_PARSER Parsing Begins..."<<endl;
   read_csv();
   run_parse();
}
///////////////////////////////////////////////////////////////////////////////
void icr_parser::make_grid(void){

   //makes the grid for evaluation of each component
   nelem = 1;
   vector<double> dimenstmp;
   for(int i=0;i<(int) grid_params.size();i++){
       dimenstmp.clear();
      if(set[i] ==0){
          for(double j=grid_params[i][0];j<grid_params[i][1];j+=grid_params[i][2]){
//cout<<j<<endl;
            dimenstmp.push_back(j);
          }
          dimens.push_back(dimenstmp);
          nelem *= dimens[i].size(); 
          activedims.push_back(i);
      } else {
          inactivedims.push_back(i);
          dimenstmp.push_back(grid_params[i][0]);
          dimens.push_back(dimenstmp);
      }
   }
//cout<<nelem<<"\t"<<ncol<<endl;
   evalpts.set_size(nelem, activedims.size());

   long ctnest=1,jctr=0,activedimctr=0;
//cout<<"dimsens.size() "<<dimens.size()<<endl;
   for(int i=0;i<(int)dimens.size();i++){
       if(set[i] == 1){
         // for(int j=0;j<nelem;j++){
         //    evalpts(i,j) = dimens[i][0];
         // }
         //Don't do anything!
       } else {
//cout<<"i:\t"<<i<<endl;
          jctr=0;
          for(int m=0;m<(int)nelem/dimens[i].size()/ctnest;m++){
             for(int k=0;k<(int)dimens[i].size();k++){
                for(int l=0;l<ctnest;l++){
//cout<<"activedimctr\t"<<activedimctr<<"\tjctr\t"<<jctr<<endl;
                    evalpts(jctr,activedimctr) = dimens[i][k];
//cout<<"jctr\t"<<jctr<<endl;
                    jctr++;
                }
             }
          }
          ctnest*=dimens[i].size();
          activedimctr++;
       }
   }
   output.set_size(nelem);
//   evalpts.print("evalpts");
//cout<<nelem<<endl;
//cout<<"AJKASKLBFAJ"<<endl;
   
}
///////////////////////////////////////////////////////////////////////////////
void icr_parser::run_parse(void){
   mat cond_mu, this_cov, cond_cov;
   mat mu1, mu2,cov11,cov12,cov21,cov22, cov22_inv;
   mat output_mat;
   rowvec a;
   long con_Nmq, con_q;
   ifstream myfile;
   myfile.open(datafile.c_str(), ifstream::binary);
   double buffer;
   if( myfile.is_open()){
     myfile.read((char *)(&buffer), sizeof(buffer));
     niter = buffer;
     myfile.read((char *)(&buffer), sizeof(buffer));
     N = buffer;
     myfile.read((char *)(&buffer), sizeof(buffer));
     ncol = buffer;
     cout<<"Niter:"<<niter<<", N:"<<N<<", ncol:"<<ncol<<endl;
     mu_star.set_size(N,ncol);
     sigma_star.set_size(ncol, ncol,N);
     njstar.set_size(N);
     make_grid();
     cout<<"Output matrix will have "<<evalpts.n_elem*ceil(niter/10)<<" elements."<<endl;
     long output_ctr=0;
     for(int i=0;i<niter;i++){
         if(i % 10 == 0 && i >0.3*niter) output_ctr++;
     }
     output_mat.set_size(evalpts.n_rows, output_ctr);
     con_q = activedims.size();
     con_Nmq = inactivedims.size();
     a.set_size(inactivedims.size());
     for(int i=0;i<inactivedims.size();i++){
       a(i) = dimens[i][0];
     }
     a.print("a:");
     permutation.set_size(ncol);
     //permutation allows me to group dimensions for division
     for(int i=0;i<(long)activedims.size();i++){
         permutation(activedims[i]) = i;
     }
     for(int i=0;i<(long)inactivedims.size();i++){
         permutation(inactivedims[i]) = i+activedims.size();
     }
     //permutation.print("permutation:");
     int iter_ctr = 0;
     while(!myfile.eof()){
         //read in each iteration
         myfile.read((char *)(&buffer), sizeof(buffer));
         iternum = buffer;
//cout<<iternum<<endl;
       for(int i=0;i<N;i++){
         for(int j=0;j<ncol;j++){
            myfile.read((char *)(&(buffer)),sizeof(buffer));
            mu_star(i,permutation[j]) = buffer;
         }
       }
//mu_star.print("mustar:");
       for(int i=0;i<N;i++){
//          mat tmp = state->sigma_star.subcube(span::all,span::all,span(i,i));
//          tmp.print("tmp");
          for(int j=0;j<ncol;j++){
            for(int k=0;k<ncol;k++){
              myfile.read((char *)(&(buffer)),sizeof(buffer));
              sigma_star(permutation[j],permutation[k],i)=buffer;
            }
           }
        }
        for(int i=0;i<N;i++){
          myfile.read((char *)(&(buffer)),sizeof(buffer));
          njstar[i]=buffer;
        }
        ndata = sum(njstar);
        myfile.read((char *)(&(buffer)),sizeof(buffer));
        alpha = buffer;
        output.zeros();
        if(con_Nmq != 0){
        for(int i=0;i<N;i++){
           this_cov = ((mat) sigma_star.subcube(span(span::all), 
                                span(span::all), span(i,i)));
           cov11 = this_cov.submat(span(0,con_q-1), span(0, con_q-1));
           cov22 = this_cov.submat(span(con_q, ncol-1), span(con_q, ncol-1));
           cov12 = this_cov.submat(span(0, con_q-1), span(con_q, ncol-1));
           cov21 = this_cov.submat(span(con_q, ncol-1), span(0, con_q-1));
           cov22_inv = pinv(cov22);
           cond_mu = mu_star.submat(span(i,i), span(0,con_q-1)) 
                   + cov12*cov22_inv*(a - mu_star.submat(span(i,i), span(con_q,ncol-1)));
           cond_cov = cov11 - cov12*cov22_inv*cov21;
/*if(iternum > 9000){
if(njstar[i] == max(njstar)){
double tmp;
tmp=0;
//cout<<iternum<<endl;
cout<<sum(njstar)<<endl;
njstar.print("njstar:");
cout<<njstar(i)<<endl;
cout<<"alpha:"<<alpha<<endl;
//cout<<njstar[i]<<endl;
           this_cov.print("this_cov:");
           mu_star.row(i).print("mu_star"); 
           //cov11.print("cov11:");
           //cov12.print("cov12:");
           //cov21.print("cov21:");
           //cov22.print("cov22:");
           //cov22_inv.print("cov22_inv:");
           // mu_star.print("mu_star:");
           cond_mu.print("cond_mu:");
           cond_cov.print("cond_cov:");
 exit(1);
 }
}*/
           output_tmp =  exp(dmvnorm(evalpts, cond_mu, cond_cov)); 
           output += output_tmp*njstar(i)/ndata;
        }} else {
        for(int i=0;i<N;i++){
           this_cov = ((mat) sigma_star.subcube(span(span::all), 
                                span(span::all), span(i,i)));
           cond_mu = mu_star.row(i);
           output_tmp =  exp(dmvnorm(evalpts, cond_mu, this_cov)); 
           output += output_tmp*njstar[i]/ndata;
        }
        }
        
        //mu_star.print("mu_star:");
        output_mat.col(iter_ctr) = output;
        iter_ctr++;
          
   }
 }
 myfile.close();
 output_mat.save("output.txt",  arma_ascii);
 evalpts.save("evalpts.txt", arma_ascii);

}
///////////////////////////////////////////////////////////////////////////////
vector<double> icr_parser::getNextLineAndSplitIntoTokens(istream& str)
{
    //vector<string>   result;
    string                line;
    getline(str,line);

    stringstream          lineStream(line);
    string                cell;
    vector<double>  output;

    while(getline(lineStream,cell,','))
    {
        output.push_back( atof(cell.c_str()));
    }
    return output;
}


///////////////////////////////////////////////////////////////////////////////

void icr_parser::read_csv(void){

  ifstream myfile(infile);
  if(myfile.is_open()){
    quantiles=getNextLineAndSplitIntoTokens(myfile);
    while(!myfile.eof()){
       grid_params.push_back(getNextLineAndSplitIntoTokens(myfile));
       if(grid_params[(int)grid_params.size()-1].size() ==1){
          set.push_back(1);
       } else {set.push_back(0);}
    }

  }
  myfile.close();
  //test_read();
}
///////////////////////////////////////////////////////////////////////////////
void icr_parser::test_read(void){
cout<<"Repeat parser input"<<endl;
 cout<<"Quantiles:"<<endl;
  for(int i=0; i<(int) quantiles.size();i++){
    cout<<quantiles[i]<<",";
  }
  cout<<endl;
 cout<<"Grid:"<<endl;
  for(int i=0; i<(int) grid_params.size();i++){
    for(int j=0; j<(int) grid_params[i].size();j++){
        cout<<grid_params[i][j]<<",";
    } 
    cout<<endl;
  }
 cout<<"Set:"<<endl;
 for(int i=0; i<(int) set.size();i++){
   cout<<set[i]<<",";
 }
 cout<<endl;
}
///////////////////////////////////////////////////////////////////////////////



icr_parser::~icr_parser(){
  //destructor
}
